﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace AcceleratorTests
{
    public static class MatrixMath
    {

        public class array2dwrapper<T>
        {
            private T[,] x;
            public array2dwrapper(int dim)
            {
                x = new T[dim, dim];
            }

            public array2dwrapper(T[,] b)
            {
                x = b;
            }

            public T this[int i1, int i2]
            {
                get
                {
                    return x[i1 - 1, i2 - 1];
                }
                set
                {
                    x[i1 - 1, i2 - 1] = value;
                }
            }


        }

        public class array1dwrapper<T>
        {
            private T[] x;
            public array1dwrapper(int dim)
            {
                x = new T[dim];
            }

            public array1dwrapper(T[] b)
            {
                x = b;
            }

            public T this[int i1]
            {
                get
                {
                    return x[i1 - 1];
                }
                set
                {
                    x[i1 - 1] = value;
                }
            }
        }

        public static void Multiply(float[,] R, float[,] A, float[,] B)
        {
            int
              aRows = A.GetLength(0),
              aCols = A.GetLength(1),
              bRows = B.GetLength(0),
              bCols = B.GetLength(1),
              rRows = R.GetLength(0),
              rCols = R.GetLength(1);
            if (aCols == bRows && rRows == aRows && rCols == bCols)
            {
                for (int r = 0; r < rRows; r++)
                    for (int c = 0; c < rCols; c++)
                    {
                        float sum = 0.0f;
                        for (int k = 0; k < aCols; k++)
                            sum += A[r, k] * B[k, c];
                        R[r, c] = sum;
                    }
            }
        }

    public static double ludcmp(float[,] a, int n, int[] index) {
    int i, imax = 0, j, k;
    float d, big, dum, sum, temp;
    float[] vv = new float[n];

    d = 1.0f;
    for (i = 0; i < n; i++) {
      big = 0.0f;
      for (j = 0; j < n; j++)
        if ((temp = Math.Abs(a[i,j])) > big) big = temp;
      if (big == 0.0) {
        Console.WriteLine ("Singular matrix in routine ludcmp");
        throw new Exception("Singular matrix in routine ludcmp");
      }
      vv[i] = 1.0f / big;
    }

    for (j = 0; j < n; j++) {
      for (i = 0; i < j; i++) {
        sum = a[i,j];
        for (k = 0; k < i; k++) sum -= a[i,k] * a[k,j];
        a[i,j] = sum;
      }
      big = 0.0f;
      for (i = j; i < n; i++) {
        sum = a[i,j];
        for (k = 0; k < j; k++) sum -= a[i,k] * a[k,j];
        a[i,j] = sum;
        if ((dum = vv[i] * Math.Abs(sum)) >= big) {
          big = dum;
          imax = i;
        }
      }
      if (j != imax) {
        for (k = 0; k < n; k++) {
          dum = a[imax,k];
          a[imax,k] = a[j,k];
          a[j,k] = dum;
        }
        d = -d;
        vv[imax] = vv[j];
      }
      index[j] = imax;
      if (a[j,j] == 0.0) a[j,j] = 1.0e-20f;
      if (j != n - 1) {
        dum = 1.0f / a[j,j];
        for (i = j + 1; i < n; i++) a[i,j] *= dum;
      }
    }

    return d;
  }




    public static float[,] inverse(float[,] a)
    {
        int i, j, N = a.GetUpperBound(0)+1;
        int[] index = new int[N];
        float[] col =  new float[N];

        float[,] y = new float[N, N];

        ludcmp(a, N, index);

        for (j = 0; j < N; j++)
        {
            for (i = 0; i < N; i++) col[i] = 0.0f;
            col[j] = 1.0f;
            lubksb(a, N, index, col);
            for (i = 0; i < N; i++) y[i,j] = col[i];
        }

        return y;
    }

   public static void lubksb(float[,] a, int n, int[] index, float[] b) {
    int i, ii = -1, ip, j;
    float sum;

    for (i = 0; i < n; i++) {
      ip = index[i];
      sum = b[ip];
      b[ip] = b[i];
      if (ii != -1) for (j = ii; j <= i - 1; j++) sum -= a[i,j] * b[j];
      else if (sum != 0.0) ii = i;
      b[i] = sum;
    }

    for (i = n - 1; i >= 0; i--) {
      sum = b[i];
      for (j = i + 1; j < n; j++) sum -= a[i,j] * b[j];
      b[i] = sum / a[i,i];
    }
  }

    }
}
